//Examples/jerasure_03.cpp

#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "jerasure.h"
#include "galois.h"

#define talloc(type, num) (type *) malloc(sizeof(type)*(num))

void usage(char *s)
{
	fprintf(stderr, "usage: jerasure_03 k w - Creates a kxk Cauchy matrix in GF(2^w). \n\n");
	fprintf(stderr, "       k must be < 2^w.  Element i,j is 1/(i+(2^w-j-1)).  (If that is\n");
	fprintf(stderr, "       1/0, then it sets it to zero).  \n");
	fprintf(stderr, "       It then tests whether that matrix is invertible.\n");
	fprintf(stderr, "       If it is invertible, then it prints out the inverse.\n");
	fprintf(stderr, "       Finally, it prints the product of the matrix and its inverse.\n");
	fprintf(stderr, "       \n");
	fprintf(stderr, "This demonstrates: galois_single_divide()\n");
	fprintf(stderr, "                   JER_Matrix::Set()\n");
	fprintf(stderr, "                   JER_Matrix::Print()\n");
	fprintf(stderr, "                   Inverse()\n");
	fprintf(stderr, "                   Prod()\n");
	if (s != NULL) fprintf(stderr, "%s\n", s);
	exit(1);
}

int main(int argc, char **argv)
{
	unsigned int k, w, i, j, n;
	JER_Matrix *jm;
	JER_Matrix *jm_prod;
	JER_Matrix *jm_inv;

	if (argc != 3) usage(NULL);
	if (sscanf(argv[1], "%d", &k) == 0 || k <= 0) usage((char *)"Bad k");
	if (sscanf(argv[2], "%d", &w) == 0 || w <= 0 || w > 31) usage((char *)"Bad w");
	if (k >= (1 << w)) usage((char *)"K too big");

	jm = new JER_Matrix(k, k, w);

	for (i = 0; i < k; i++) {
		for (j = 0; j < k; j++) {
			n = i ^ ((1 << w)-1-j);
			jm->Set(i,j,(n == 0) ? 0: galois_single_divide(1, n, w));
		}
	}

	printf("The Cauchy Matrix:\n");
	jm->Print();
	
	jm_inv = Inverse(jm);
	if(jm_inv != NULL){
		printf("\nInvertible: Yes\n");

		printf("\nInverse:\n");
		jm_inv->Print();

		printf("\nInverse times matrix (should be identity):\n");
		jm_prod = Prod(jm_inv,jm);
		
		if(jm_prod != NULL){
			jm_prod->Print();	
			delete jm_prod;
		}else{
			printf("\nError\n");
		}
	}else{
		printf("\nInvertible: No\n"); 
	}

	delete jm_inv;
	delete jm;

	return 0;
}

